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ABSTRACT 

WU,  W.;  SANCHEZ,  A.,  and  ZHANG,  M„  2011.  An  Implicit  2-D  Shallow  Water  Flow  Model  on  Unstructured 
Quadtree  Rectangular  Mesh.  In:  Roberts,  T.M.,  Rosati,  J.D.,  and  Wang,  P.  (eds.),  Proceedings,  Symposium  to  Honor 
Dr.  Nicholas  C.  Kraus,  Journal  of  Coastal  Research,  Special  Issue,  No.  59,  pp.  15-26.  West  Palm  Beach  (Florida), 
ISSN  0749-0208. 

An  implicit  finite  volume  scheme  is  developed  to  solve  the  depth-averaged  2-D  shallow  water  flow  equations.  The 
computational  mesh  consists  of  rectangular  cells,  with  quadtree  technology  incorporated  to  locally  refine  the  mesh 
around  structures  of  interest  or  where  the  topography  and/or  flow  properties  change  sharply.  The  grid  nodes  are 
numbered  by  means  of  an  unstructured  index  system  for  more  flexibility.  The  governing  equations  are  solved  using 
the  SIMPLEC  algorithm  on  non-staggered  grid  to  handle  the  coupling  of  water  level  and  velocity.  In  this  non- 
staggered  system,  primary  variables  u-,  v-velocity,  and  water  level  are  stored  on  the  same  set  of  grid  points,  and 
fluxes  at  cell  faces  are  determined  using  the  Rhie  and  Chow’s  momentum  interpolation  method  to  avoid  spurious 
checkerboard  oscillations.  The  discretized  algebraic  equations  are  solved  iteratively  using  the  GMRES  method.  The 
model  has  been  tested  against  measurement  data  for  steady  flow  around  a  spur-dyke  in  a  laboratory  flume  and  tidal 
flows  in  Gironde  Estuary,  France  and  Grays  Harbor,  USA.  The  model  reasonably  well  reproduces  the  temporal  and 
spatial  variations  of  water  level  and  current  speed  observed  in  the  measurements.  The  laboratory  test  has 
demonstrated  that  the  quadtree  mesh  is  cost-effective,  while  the  two  field  cases  have  shown  that  the  model  is  very 
stable  and  handles  wetting  and  drying  efficiently. 

ADDITIONAL  INDEX  WORDS:  Shallow  water  flow  equations,  two-dimensional,  finite  volume,  numerical  model, 
quadtree  rectangular  mesh. 


INTRODUCTION 

Because  of  their  nonlinearity  and  irregular  domains,  most  of 
the  real-life  problems  of  surface  water  flows  in  rivers,  lakes, 
estuaries  and  coastal  water  bodies  have  to  be  solved 
numerically.  The  numerical  methods  widely  used  include  finite 
difference  method  (FDM)  (Fletcher,  1991),  finite  element 
method  (FEM)  (Chung,  1978;  Zienkiewicz  and  Taylor,  2000), 
finite  volume  method  (FVM)  (Patankar,  1980;  Ferziger  and 
Peric,  1995;  Wu,  2007),  etc.  The  algebraic  equations  resulting 
from  the  classic,  structured  FDM  and  FVM  usually  have  banded 
coefficient  matrices  that  can  be  handled  efficiently,  whereas  the 
algebraic  equations  from  the  unstructured  FEM  usually  have 
sparse  coefficient  matrices  that  require  relatively  tedious  effort 
for  solution.  However,  the  classic  FDM  and  FVM  usually  adopt 
regular  meshes  and  encounter  difficulties  in  conforming  to  the 
irregular  domains  of  surface  water  flows,  while  the  FEM 
typically  uses  irregular  meshes  that  can  conveniently  handle 
such  irregular  domains.  Therefore,  it  has  been  a  trend  in  recent 
decades  to  develop  the  FDM  and  FVM  based  on  structured  or 
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unstructured  irregular  meshes  with  quadrilaterals,  triangles  and 
polygons,  which  have  the  grid  flexibility  of  the  FEM  and  the 
merits  of  the  classic  FDM  and  FVM. 

Enhancement  of  accuracy  of  numerical  solution  is  one  of  the 
main  concerns  in  computational  simulation  of  free  surface 
flows.  One  may  simply  use  high-order  accurate  schemes  to 
discretize  the  differential  governing  equations  for  this  purpose, 
but  a  high-order  scheme  involves  more  computational  nodes,  is 
more  complicated  and  requires  more  computational  time. 
Numerical  schemes  higher  than  third  or  fourth  order  have  been 
rarely  used  in  simulation  of  surface  water  flows.  Another 
approach  to  enhancing  accuracy  is  through  refinement  of  mesh. 
Numerical  discretization  based  on  a  finer  mesh  usually  has  less 
truncation  errors  and  thus  yields  more  accurate  solution. 
However,  globally  refining  a  large,  complex  computational 
mesh  is  often  costly,  and  only  locally  refining  near  boundaries 
and  high-gradient  regions  is  sometimes  needed.  For  local 
refinement,  an  unstructured  triangular  mesh  is  a  good  choice, 
whereas  a  structured  rectangular  mesh  is  inconvenient  because  it 
refines  the  mesh  in  regions  which  are  unnecessary.  A  structured 
quadrilateral  (curvilinear)  grid  can  serve  this  purpose  by 
stretching  or  shrinking  the  mesh  lines,  but  it  is  less  flexible  for 
very  large,  complex  domains  than  the  triangular  mesh.  On  the 
other  hand,  the  rectangular  or  quadrilateral  mesh  is  more 
convenient  than  the  triangular  mesh  for  establishing  high-order 
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(e.g.,  second  and  third)  schemes  or  for  discretizing  second  or 
higher  order  spatial  derivatives.  Therefore,  quadtree 
(telescoping)  rectangular  or  quadrilateral  mesh  has  been  recently 
used  for  local  refinement  of  computational  mesh  in 
computational  fluid  dynamics  (CFD)  (Muzaferija,  1994; 
Ferziger  and  Peric,  1995;  Greaves,  2004;  Liang  et  al.,  2008; 
Nabi,  2008).  On  the  quadtree  mesh,  a  coarse  cell  is  split  into 
four  child  cells,  and  as  many  levels  of  refinement  as  necessary 
can  be  used.  The  quadtree  mesh  can  be  arranged  as  block- 
structured  or  completely  unstructured.  It  is  flexible  in  mesh 
generation,  while  keeping  most  advantages  of  the  rectangular  or 
quadrilateral  mesh. 

For  incompressible  flows  such  as  surface  water  flows,  the 
governing  equations  are  the  Navier-Stokes  equations,  in  which 
the  velocity  is  linked  to  the  pressure  gradient  by  the  momentum 
equations  but  not  the  continuity  equation.  The  continuity 
equation  is  just  an  additional  constraint  on  the  velocity  field 
without  directly  linking  to  the  pressure.  Because  of  such  a  weak 
linkage,  the  convergence  and  stability  of  a  numerical  solution  of 
the  Navier-Stokes  equations  depends  largely  on  the  coupling 
between  the  pressure  and  velocity  fields.  Storing  the  variables  at 
the  geometric  center  of  the  control  volume  in  combination  with 
linear  interpolation  for  intemodal  variation  usually  leads  to  non¬ 
physical  node-to-node  (checkerboard)  oscillations.  One 
approach  for  eliminating  such  oscillations  is  to  use  the  staggered 
grid,  as  adopted  in  Harlow  and  Welch’s  (1965)  MAC  (Marker 
and  Cell)  method,  Chorin’s  (1968)  projection  method,  and 
Patankar  and  Spalding’s  (1972)  SIMPLE  (Semi-Implicit  Method 
for  Pressure-Linked  Equations)  algorithm.  The  other  approach  is 
to  use  the  momentum  interpolation  technique  proposed  by  Rhie 
and  Chow  (1983)  based  on  the  non-staggered  grid.  For  the 
depth-averaged  2-D  simplification  of  the  Navier-Stokes 
equations  in  the  case  of  shallow  water  flows,  the  linkage 
between  the  flow  velocity  and  pressure  (water  level)  is  improved 
in  the  continuity  equation,  but  the  ware-surface-gradient  terms 
remain  in  the  momentum  equations  and  the  issue  of  coupling 
velocity  and  water  level  still  exists  somehow.  In  recent  years, 
the  staggered  grid  approaches  and  the  Rhic  and  Chow’s  (1983) 
momentum  interpolation  method  on  non-staggered  grid  have 
been  applied  to  the  depth-averaged  simulation  of  surface  water 
flows  (Wenka,  1992;  Lu  and  Zhang,  1993;  Ye  and 
McCorquodale,  1997;  Minh  Due,  1998;  Kim  et  al,  2003;  Wu, 
2004).  In  addition,  several  attempts  {e.g.,  Lc  Roux  et  al.,  1998) 
have  been  made  in  the  FEM  to  eliminate  the  non-physical 
oscillations. 

This  paper  presents  recent  advances  in  the  Coastal  Modeling 
System  (CMS)  for  nearshore  circulation  developed  under  the 
Coastal  Inlets  Research  Program  (CIRP)  of  U.S.  Army  Corps  of 
Engineers  (Militello  et  al.,  2004;  Buttolph  et  al.,  2006).  The 
existing  CMS  circulation  model  (called  CMS-Flow)  solves  the 
2-D  shallow  water  equations  using  an  explicit  FVM  on 
structured  rectangular  mesh.  It  computes  the  unsteady  water 
level  and  current  velocity  fields  by  solving  the  depth-averaged 
2-D  shallow  water  flow  equations  on  a  non-uniform  Cartesian 
grid  with  an  explicit  finite  volume  scheme.  It  can  simulate  tide, 
wind  and  wave  driven  currents.  It  is  two-way  coupled  with  the 
spectral  wave  transformation  model  called  CMS-Wave,  which 
solves  the  steady-state  wave-action  balance  equation  on  a  non- 
uniform  Cartesian  grid  with  a  finite  difference  scheme  and 


considers  wind  wave  generation  and  growth,  diffraction, 
reflection,  dissipation  due  to  bottom  friction,  white  capping  and 
breaking,  wave-wave  and  wave-current  interactions,  wave 
runup,  wave  setup,  and  wave  transmission  through  structures 
(Lin  et  al.,  2008).  To  improve  the  CMS-Flow  model’s 
computational  efficiency,  an  implicit  depth-averaged  2-D 
scheme  has  been  designed  and  tested  to  solve  the  shallow  water 
flow  equations  based  on  a  quadtree  rectangular  mesh.  This 
model  solves  the  2-D  shallow  water  equations  for  arbitrary 
combinations  of  currents  and  waves,  considering  the  effect  of 
phase-averaged  wave  radiation  stresses  on  the  current.  The 
governing  equations  are  solved  within  the  FVM  with  the 
SIMPLEC  (SIMPLE  Consistent;  van  Doormaal  and  Raithby, 
1984)  algorithm  on  non-staggered  grid  to  handle  the  coupling  of 
water  level  and  velocity.  Due  to  limit  of  paper  length,  the 
numerical  solution  methodologies  used  in  the  new  CMS  flow 
model  are  the  main  focus  here,  while  the  wave/current 
interaction,  bottom  friction,  cross-shore  boundary  condition, 
sediment  and  salinity  transport,  and  many  other  related  issues 
will  be  addressed  on  other  occasions. 


GOVERNING  EQUATIONS  AND  BOUNDARY 
CONDITIONS 


The  depth-averaged  shallow  water  equations  in  the  Cartesian 
coordinate  system  are  written  as 


8h  8{hu )  +  d(hv) 
8t  dx  dy 


8{hu )  d(huu  )  8(huv ) 
8t  dx  dy 


d?l  5  8  f.,u8u 


-gh — -  H - v.h —  H - v.h — 

dx  dx  v  dx )  dyy  dy 


+—(TS,+Tw,-Tb1)+fchv 

P 

8{hv )  d(huv)  d(hvv) 
dt  dx  dy 


Tsy+Tv-Tby)-fchU 


CD 


(2) 


0) 


where  t  is  the  time;  x  and  y  are  the  horizontal  Cartesian 
coordinates;  h  is  the  total  flow  depth;  u  and  v  are  the  depth- 
averaged  flow  velocities  in  x-  and  y-directions;  it  is  the  water 
surface  elevation  above  the  reference  sea  level;  g  is  the 
gravitational  acceleration;  p  is  the  density  of  flow;  v,  is  the  eddy 
viscosity  due  to  turbulence;  rbx  and  zby  are  the  bed  shear  stresses 
in  x-  and  y-directions  that  are  determined  by  =  pCfu\lu2  +  v2 

and  r by=pCfv4urW  J  cf=gn2/h'n  .  in  which  n  is  the 

Manning’s  roughness  coefficient;  Tgx  and  zyv,  are  the  wave 
radiation  stresses  in  x-  and  y-directions;  twx  and  Twy  are  the  wind 
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driving  forces;  and  f.  is  the  Coriolis  force  coefficient.  Because 
the  present  paper  focuses  mainly  on  the  numerical  solution 
methods  of  the  established  flow  model,  detennination  of  the 
wind  driving  force  and  the  wave  radiation  stresses  refers  to 
Buttolph  et  al.  (2006)  and  Lin  et  al.  (2008). 

Several  turbulence  models,  including  the  depth-averaged 
parabolic  eddy  viscosity  model  and  the  modified  mixing  length 
model,  have  been  implemented  in  the  developed  model  to 
detennine  the  eddy  viscosity  v,.  In  the  depth-averaged  parabolic 
model,  the  eddy  viscosity  is  calculated  by  vt=au*h,  in  which  u*  is 
the  bed  shear  velocity,  «»=[cXtt2+v2)]1/2,  and  a  is  an  empirical 
coefficient  between  0.3-1 .0.  It  is  noted  that  if  a  is  set  to  0.578cy, 
then  the  depth-averaged  parabolic  model  reduces  to  the  Falconer 
(1980)  equation  for  eddy  viscosity  previously  used  in  the  CMS- 
Flow  model. 

The  depth-averaged  parabolic  eddy  viscosity  model  is 
applicable  in  the  region  of  main  flow  but  does  not  account  for 
the  influence  of  the  horizontal  gradient  of  velocity. 
Improvement  can  be  achieved  through  combination  of  the  depth- 
averaged  parabolic  eddy  viscosity  model  and  the  mixing-length 
model  (Wu,  2007): 


v,  = 


(a0ujif  +(/j  l^l)2 


(4) 


where  |s|  =  ^2(du/dx)2  +2(5v/5y)2  +  (du/dy  +  dv/dx)~  J  >  ao  is 

an  empirical  coefficient,  set  as  rc/6,  with  k  being  the  von  Kannan 
constant;  lh  is  the  horizontal  mixing  length,  determined  by 
/  =  Km\n(cji,  y)  ,  with  y  being  the  distance  to  the  nearest  wall 

and  c,„  an  empirical  coefficient  between  0.3-1. 2.  The  coefficient 
cm  is  often  given  larger  values  in  field  cases  and  smaller  values 
in  laboratory  cases,  ao  is  smaller  than  a,  because  Equation  (4) 
takes  into  account  the  effect  of  horizontal  velocity  gradients 
through  the  second  term  on  its  right-hand  side. 

For  a  well-defined  problem  governed  by  Equations  ( 1  )-(3), 
the  flow  discharge  or  velocity  is  needed  at  inflow  boundaries, 
while  the  water  level  is  usually  given  at  outflow  boundaries  for  a 
subcritical  flow  or  at  inflow  boundaries  for  a  supercritical  flow. 
Near  rigid  wall  boundaries,  such  as  beaches  and  islands,  the 
wall-function  approach  is  employed.  By  applying  the  log-law  of 
velocity,  the  resultant  wall  shear  stress,  f  ,  is  related  to  the  flow 

velocity,  y  ,  at  the  center,  P,  of  the  control  volume  close  to  the 
wall  by  the  following  relation: 


,=-Wv 


(5) 


where  X  is  a  coefficient  detennined  as  x  =  putKj\a(Ey+p^  with 

y*P  =  pu,y pj //>  'n  which  fi  is  the  dynamic  viscosity,  yp  is  the 

distance  from  cell  center  P  to  the  wall,  and  £  is  a  coefficient 
related  to  wall  roughness  (Wu,  2007).  Since  X  is  related  to  u*, 
iteration  is  needed  to  solve  Equation  (5). 


QUADTREE  GRID  AND  DATA  STRUCTURE 

Because  of  the  complexity  of  computational  domain  near 
coastal  inlets,  a  simple  structured  rectangular  mesh  requires  a 
large  number  of  cells  to  resolve  the  detailed  flow  pattern  near 
inlets,  navigation  channels  and  in-stream  structures.  To  optimize 
the  use  of  computational  resources,  we  use  the  multiple-level 
quadtree  rectangular  mesh  with  local  refinement.  In  this  mesh, 
various  levels  of  fine  cells  are  placed  close  to  the  navigation 
channels  and  in-stream  structures  where  the  flow  gradients  are 
high,  while  coarse  grids  are  used  in  the  low-gradient  regions. 
For  simplifying  the  mesh,  a  cell  is  refined  by  splitting  into  four 
equal  child  cells.  Corresponding  to  this  refining,  any  cell  has  one 
or  two  faces  on  each  of  its  south,  north,  west,  and  east  sides.  For 
further  simplification,  we  eliminate  those  isolated  single  refined 
or  coarse  cells.  This  means  that  a  cell  should  be  refined  if  all  of 
its  adjacent  cells  on  either  x  or y  direction  are  refined,  and  on  the 
other  hand,  a  cell  should  not  be  refined  if  all  of  its  adjacent  cells 
are  not  refined.  Through  this  handling,  each  cell  has  only  four  to 
six  faces  even  though  its  each  side  may  have  one  or  two  faces, 
as  shown  in  Figures  1  and  2,  so  that  the  computational  mesh  will 
be  less  complicated. 

The  data  structure  for  the  quadtree  mesh  can  be  managed  in 
several  ways:  block-structured,  hierarchical  tree,  and  fully 
unstructured.  The  block-structured  approach  divides  the  domain 
into  multiple  blocks,  each  of  which  is  treated  as  structured. 
Interfaces  between  blocks  are  specially  handled  to  ensure  mass 
and  momentum  balance  through  them.  The  tree  data  structure 
uses  parent  and  child  relations  and  requires  tree  traverse  to 
determine  the  mesh  connectivity.  In  the  fully  unstructured 
approach,  all  cells  are  numbered  in  a  one-dimensional  sequence, 
and  pointers  are  used  to  detennine  the  connectivity  of 
neighboring  cells  for  each  cell.  Among  the  three  approaches,  the 
fully  unstructured  approach  is  simpler  and  thus  is  used  in  this 
study. 


Figure  1 .  Example  of  quadtree  mesh  system. 
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Figure  2.  Control  volume  in  a  quadtree  mesh. 


As  mentioned  in  the  Introduction,  another  issue  in  simulation 
of  incompressible  flow  is  the  location  of  primary  variables: 
velocity  and  pressure  (water  level),  related  to  elimination  of  the 
non-physical  node-to-node  (checkerboard)  oscillations.  One  can 
use  either  staggered  or  non-staggered  grid.  On  a  staggered  grid, 
the  pressure  is  located  at  the  center  of  cells  and  the  u-  and  v- 
velocities  are  on  the  faces  of  cells  (Harlow  and  Welsh,  1965; 
Patankar,  1980).  On  the  non-staggered  grid,  all  the  primary 
variables  are  located  at  the  center  of  cells.  The  staggered  grid 
can  more  conveniently  eliminate  the  checkerboard  oscillations 
than  the  non-staggered  grid,  but  the  non-staggered  grid  results  in 
a  simpler  computer  code  and  can  minimize  the  number  of 
coefficients  that  must  be  computed  and  stored  because  many  of 
the  tenns  in  each  of  the  equations  are  essentially  identical.  In 
particular,  the  staggered  grid  is  more  complicated  in  handling 
the  interface  between  coarse  and  fine  cells  where  five-  or  six- 
face  control  volumes  are  used.  Therefore,  the  non-staggered  grid 
approach  is  adopted  here,  with  the  Rhic  and  Chow’s  momentum 
interpolation  technique  used  to  eliminate  the  checkerboard 
oscillations. 


NUMERICAL  DISCRETIZATION 

Integrating  the  continuity  equation  (1)  over  the  control 
volume  shown  in  Figure  2,  applying  Green’s  theorem  and 
discretizing  the  temporal  derivative  by  the  backward  difference 
scheme,  one  can  derive  the  following  equation: 

hn+l-h"  A,.  \  n+1  ,  y,+l. 

A AP+X{hu)ek  A yet  4>U 

At  k= i  k= i  (6) 

mn  ms 

+z(/n;)r‘  a**  -  z(/nC  ax* = o 

t=i  *=i 


where  At  is  the  time  step  length;  A AP  is  the  area  of  the  control 
volume  (cell)  at  node  P;  Ax  and  Ay  are  the  lengths  of  the  cell 
faces  in  either  x  or  y  direction;  the  subscripts  w,  e,  s  and  n  denote 
the  west  (negative  x),  east  (positive  x),  south  (negative  y)  and 
north  (positive  y)  sides  of  the  control  volume;  the  subscript  k  is 


the  index  of  the  faces  at  each  side,  with  a  value  of  1  or  2;  and 
mw,  me,  ms  and  m„  are  the  numbers  of  cell  faces  at  west,  east, 
south  and  north  sides  of  the  cell.  For  the  control  volume  shown 
in  Figure  2,  mw=  1,  me=  2,  ms=  1  and  m„= 2.  For  simplicity,  mw,  me, 
ms,  m„,  and  the  superscript  n+ 1  will  be  omitted  in  the  following 
notations. 

Defining  the  fluxes  at  cell  faces  as  Fek=(hu)ekAyek, 
Fwk=Qiu)wkAywh  Fnk=(hv)nk Axnk  and  Fsk={hv)sk Axsk,  one  can 
rewrite  Equation  (6)  as 

}f^AAp  +  ^Fek-YjFwt+^Fnk-YjFsk=  0  (7) 

At  k  k  k  k 


Integration  of  the  x-momentum  equation  over  the  control 
volume  shown  in  Figure  2  leads  to 


~(hu)"p 

At 


AAP  +  Ayek  ~  YXhuu)Wk  Ay* 


+  2>“VL  A*nk  -  Z(/"'v)1*  Axsk 


-ghp 


'Zjflektyek 


du 


A Av, 


du 


dx 


(8) 


where  5„  includes  all  the  remaining  terms.  The  convection  tenns 
can  be  discretized  using  several  numerical  schemes  with 
upwinding  capability,  such  as  the  hybrid  upwind/central  scheme 
(Spalding,  1972),  exponential  scheme  (Spalding,  1972)  and 
HLPA  scheme  (Zhu,  1991),  while  the  diffusion  terms  are 
discretized  using  the  central  difference  scheme.  The  hybrid 
scheme  uses  the  first-order  upwind  or  second-order  central 
difference  scheme  depending  on  whether  the  Peclet  number  is 
larger  than  2  or  not.  The  exponential  scheme  uses  the  analytical 
solution  of  the  linearized  1-D  convection-diffusion  equation  to 
approximate  the  profile  of  the  sought  quantity  between  two 
neighboring  cell  centers;  its  accuracy  is  between  first  and  second 
orders  in  general.  The  HLPA  scheme  uses  a  combination  of 
linear  and  parabolic  approximations  for  such  profile,  and  it  has 
second-order  accuracy.  Normally  the  HLPA  scheme  needs 
slightly  more  computational  time  and  special  treatments  near  a 
boundary  and  during  the  iteration  because  it  uses  more 
computational  nodes.  The  details  can  be  found  in  Wu  (2007).  In 
cases  where  higher  accuracy  is  concerned,  the  HLPA  scheme  is 
a  better  option.  However,  for  most  of  application  cases,  the 
exponential  scheme  is  recommended  because  it  uses  about  10- 
20%  less  computational  time  than  the  HLPA  scheme  and  has 
enough  accuracy. 

After  the  above  manipulations,  the  x-momentum  equation  (8) 
reads 
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(9) 


-OpUp  -  ghp\Y^tikAyek~Y,riwlAy^Y  S„AAp 

where  am,  ctEk,  ash  am  and  flp  are  coefficients.  Similarly,  one 
can  discretize  the  y-momentum  equation  as 


m;-Mp 


At 

^  aWkVWk  +  aEkVEk  +  ^  aNkVNk  +  aSkVSk 
k  k  k  k 

-aPvP  ~  ghP  f  ^  rjnkAxnk  -  £  r/skAxsk  )  +  SVAAP 


AAP 


(10) 


COUPLING  OF  VELOCITY  AND  WATER  LEVEL 

From  the  discretized  momentum  equation  (9),  one  can  derive 
the  following  equation  for  velocity  Up+I  '■ 


relation  between  the  water  level  and  velocity  corrections  from 
Equation  (12): 


(13) 


where  if  is  the  water  level  correction  rj'  =  (7-77*  ■  In  the 


SIMPLEC  algorithm. 


and 


The  relation  of  water  level  and 


velocity  corrections  for  the  SIMPLE  algorithm  is  similar  to 
Equation  (13),  with  if  and  if  replaced  by  D\.k  and  D\k  ■ 
Similarly,  one  can  have  the  v-equation  and  the  corresponding 
correction  equation: 


Vp  =  ocv  h;p -Y.DfX.k +(1  -av)v°p  (14) 


ZAk'*  -ZAA,'* 


(15) 


<  =Y( z«x+1 + 1  -  z Di^k + z  aa*  o- 1 ) 

aP\i  J  k  t 

where  a"  denotes  the  coefficients  for  ((-equation, 
D'ek  =  ghPAyet  /  aup  and  Dlwk  =  ghpAywk  /  a“p  ■  Note  that  the  first 
summation  in  Equation  (11)  is  applied  with  the  index,  /, 
sweeping  over  all  the  neighboring  cells  of  cell  P,  and  the 
temporal  tenn  on  the  left-hand  side  of  Equation  (9)  is  arranged 
in  Stl  and  a“  in  Equation  (1 1)  as  suggested  by  Patankar  (1980). 

Equation  (11)  is  used  to  compute  the  velocity  for  an  assumed 
water  level  field  in  an  iterative  manner.  Application  of  under¬ 
relaxation  (Majumdar,  1988)  leads  to 

k=<*u[h'p- (12) 

where  rf  is  the  guessed  water  level,  «*  is  the  approximate 
solution  of  ((-velocity,  u°p  is  the  ((-velocity  in  the  previous 
iteration  step,  H\p  denotes  the  first  tenn  on  the  right-hand  side 
of  Equation  (11)  and  au  is  the  relaxation  factor  that  is  set  as 
about  0.8  in  this  study. 

There  are  several  approaches  to  conduct  the  iteration  and 
couple  the  flow  velocity  and  water  level,  including  the  SIMPLE 
(Patankar  and  Spalding,  1972),  SIMPLEC  (van  Doonnaal  and 
Raithby,  1984),  and  IDEAL  (Sun  et  al.,  2009)  algorithms.  The 
SIMPLEC  algorithm  is  used  here  since  it  has  been  tested  well 
for  shallow  water  flow  modeling  (Wu,  2007).  One  can  derive  the 


where  H'lp  +S^ja).  >  D;k  =  E^j\\ -a^a] 

and  D2sk  =Dlj[\-a^la'p 

D2,  =  sh„Ax  ,  I  al-  Here,  av  is  the  relaxation  factor  for  the  v- 

sk  o  P  sk  P  ’ 

equation. 

It  seems  that  the  water  surface  elevation  could  be  calculated 
from  the  discretized  continuity  equation  (7),  but  in  fact  this 
might  lead  to  spurious  numerical  oscillations  for  the  collocated 
arrangement,  as  explained  by  Patankar  (1980)  and  others.  In 
order  to  avoid  the  checkerboard  splitting  for  the  collocated 
arrangement,  we  adopt  the  momentum  interpolation  technique 
proposed  by  Rhic  and  Chow  (1983)  to  evaluate  the  variable 
values  at  cell  faces  from  the  quantities  at  cell  centers.  For 
example,  the  ((-velocity  at  w-face  and  the  v-velocity  at  s-face  are 
determined  as 

Uwk  =  au  [( 1  —  fx.p  )  ,pwk  +  fx.pH^P  J 

+«„[(! -L,p)/apm  +  fx,P  ! ap^gh^Ay^tJm  ~Vp)  <16> 

+  ( ■ 1  -  «„  )  [( ■ 1  -  fx ,p  )  Kk  +fx.pU°P  ] 

Vsk  =  «v  [(l  -  fy.p  )H2 ,PSk  +  fy.pK.P  ] 

+«v  [(l  -  fy.p  )  !  aPSk  +  fy.p  1  a'p  ]  SKkEXsk  {fSk-Vp)  (  1  7) 

+  (1-«v)[(1-/v,p)4+/v,Q'p] 


•  with  Dlk  =  ghPAx„k  /  avP  and 
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and  the  corresponding  velocity  corrections  as 


Q„k  (?Jwk  }7p) 

(18) 

Qsk  (  Vsk  —  7lp  ) 

(19) 

outer  iteration  loops.  The  inner  iteration  is  designed  for  solving 
each  of  the  discretized  momentum  equations  (12)  and  (14)  and 
the  water-level-correction  equation  (22)  with  an  iteration  solver 
as  discussed  in  the  next  paragraphs.  The  outer  iteration  loop 
visits  the  u,  v  and  //'equations  in  the  following  sequence  in  each 
time  step  as  required  by  the  SIMPLEC  algorithm: 


where 

Qlk  =  [(!  -  fx,P )  I  <h  +  fx.P  1  aPK  ]  SKA >’„*  / 

i  -  «„  (i  -  fx,p )  fz 1 a“  /a“p 1  _  a'A,P  ( Z  ai  I a“p  1 

V  I  Jwk  V  /  Jp 

and 

Qsk  =  [(!  -  fy.p  )  /  am  +  f y,p  1  an  ]  sK ^sk  / 

i-o,a-/^)[Zfl,7a;j  -a,/v.^Ea/Vapj  ’ 

in  which  and  yj  are  the  weighting  factors  used  to 

interpolate  the  values  of  a  variable  at  cell  faces  w  and  .v  from  the 
values  at  two  adjoining  cell  centers  P  and  W  or  P  and  S, 
respectively;  a“pw  and  a'PS  stand  for  a“p  and  ap  when  applying 
Equations  (9)  and  (10)  on  the  cells  centered  by  W  and  S, 
respectively. 

With  the  definition  of  fluxes  at  cell  faces  and  Equations  (18) 
and  (19),  one  can  derive  the  flux  corrections  at  w  and  s  faces: 

Kk=Kk+amWm-V'p)  (2°) 

Fsk=F'k+<Wsk->l’p)  (21) 

where  <  =  auQ'wkhwkAywk  ,  ansk=arQ;khsk AJU  and  K  and  K 
are  the  fluxes  at  faces  w  and  S'  in  tenns  of  the  velocities  «*  and 
v*  evaluated  using  the  Rhie  and  Chow’s  momentum 
interpolation  method. 

Inserting  Equations  (20)  and  (21)  into  (7)  leads  to  the 
following  equation  for  water  level  correction: 

apr]p  =  2 ]  aEkrlEk  +  2j  aWkTlwk  +  2 ]  aNkrlNk  +  Z  aSk*lsc  + 
k  k  k  k 

where  anp  =^<*7  +  Mp/At  .  and 

/ 

Sn  =  ~{n; -n'p)&Ap/&t-{ 2%  -  SC  +Z^  -2X ]  • 

V  k  k  k  k  ) 


(1 )  Guess  the  water  level  field  rt f ; 

(2)  Solve  the  momentum  equations  (12)  and  (14)  to  obtain 
up  and  vp  5 

(3)  Use  the  Rhie  and  Chow’s  momentum  interpolation  to 
determine  the  velocities  and  fluxes  at  cell  faces; 

(4)  Calculate  jj'  using  Equation  (22); 

(5)  Correct  //by  rj  =  rf  +rj' ,  and  update  uP  and  vP  using 
Equations  (13)  and  (15)  and  fluxes  using  Equations 
(20)  and  (21); 

(6)  Treat  the  corrected  water  level,  //,  as  a  new  guess,  rf , 
and  repeat  the  procedure  from  steps  2  to  6  until  a 
converged  solution  is  obtained. 

Because  the  cells  on  the  quadtree  mesh  are  numbered  in  an 
unstructured  fonn,  the  coefficient  matrices  in  the  discretized 
momentum  equations  and  the  water-level-correction  equation 
are  sparse.  Selection  of  iterative  solution  solvers  for  these 
algebraic  equations  is  the  key  issue  concerning  the  overall 
performance  of  the  model.  After  a  lot  of  testing,  we  have  chosen 
a  variant  of  the  GMRES  (generalized  minimum  residual) 
method  (Saad,  1993)  to  solve  the  algebraic  equations.  The 
original  GMRES  method  (Saad  and  Schultz,  1986)  uses  the 
Amoldi  process  to  reduce  the  coefficient  matrix  to  the 
Hcssenburg  fonn  and  minimizes  at  every  step  the  norm  of  the 
residual  vector  over  a  Krylov  subspace.  The  variant  of  the 
GMRES  method  recommended  by  Saad  (1993)  allows  changes 
in  the  preconditioning  at  every  iteration  step.  We  use  the  ILUT 
(Incomplete  LU  Factorization;  Saad,  1994)  as  the  preconditioner 
to  speed  up  the  convergence.  The  details  of  this  iteration  solver 
can  be  found  in  the  above  references. 

Because  the  solution  of  each  of  the  variables  u,  v  and  rf  needs 
the  updated  intermediate  values  of  other  variables,  it  is  not 
necessary  to  reach  a  convergent  solution  in  each  inner  iteration. 
However,  an  approximately  convergent  solution  for  the  rf 
equation  (22)  will  ensure  mass  conservation  at  each  iteration 
step.  It  is  a  good  practice  to  set  more  iteration  steps  when 
solving  Equation  (22).  The  optimal  inner  and  outer  iteration  step 
numbers  are  related  to  the  iterative  solution  method,  mesh  size 
and  nature  of  the  problem.  For  most  cases  with  mesh  of  less  than 
100,000  cells,  if  the  GMRES  method  is  used,  5,  5  and  10  inner 
iteration  steps  have  been  found  to  be  good  choice  for  Equations 
(12),  (14)  and  (22),  respectively,  and  20  outer  iteration  steps  are 
usually  enough. 

WETTING  AND  DRYING  TECHNIQUES 


SOLUTION  OF  DISCRETIZED  EQUATIONS 

Because  the  dicretized  equations  are  non-linear,  they  need  to 
be  solved  iteratively.  The  iteration  process  consists  of  inner  and 


In  the  numerical  simulation  of  surface  water  flows  with 
sloped  beaches,  barriers,  spits  and  islands,  the  water  edges 
change  with  time,  with  part  of  the  nodes  being  possibly  wet  or 
dry.  In  the  present  model,  a  threshold  flow  depth  (a  small  value 
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such  as  0.02  m  in  field  cases)  is  used  to  judge  drying  and 
wetting.  If  the  flow  depth  on  a  node  is  larger  than  the  threshold 
value,  this  node  is  considered  to  be  wet,  and  if  the  flow  depth  is 
lower  than  the  threshold  value,  this  node  is  dry.  Because  a  fully 
implicit  solver  is  used  in  the  present  model,  all  the  wet  and  dry 
nodes  participate  in  the  solution.  Dry  nodes  are  assigned  a  zero 
velocity.  On  the  water  edges  between  the  dry  and  wet  nodes,  the 
wall-function  approach  is  applied. 

MODEL  TESTING 


flows  around  the  spur-dyke.  The  velocity  profiles  along  these 
cross-sections  among  the  three  simulations  have  only  very  little 
difference.  This  indicates  that  the  quadtree  mesh  technology  is 
cost-effective.  It  has  been  observed  that  the  sharp  velocity 
gradients  at  the  edge  between  the  main  and  recirculation  flow 
zones  are  under-predicted  in  all  three  runs.  This  is  due  to  the 
zero-order  turbulence  model  used.  In  order  to  capture  such 
velocity  gradients,  a  higher-order  turbulence  closure  such  as  the 
k-s  turbulence  model  is  preferred  (Wu  et  al.,  2004). 


The  developed  model  has  been  recently  enhanced  to  take  into 
account  the  wave/current  interactions  and  to  simulate  sediment 
transport  and  morphology  changes.  It  has  been  tested 
extensively  in  many  laboratory  and  field  cases.  Because  the 
main  focus  of  this  paper  is  the  flow  solver,  only  three  test  cases 
are  presented  here.  The  first  case  is  a  steady  flow  around  a  spur- 
dyke  in  a  laboratory  flume.  This  case  is  designed  to  verify  the 
developed  numerical  algorithms  based  on  quadtree  meshes.  The 
second  and  third  cases  are  the  tidal  flows  in  estuaries,  through 
which  the  stability,  efficiency  and  reliability  of  the  model  for 
unsteady  flows  are  quantitatively  validated. 

Case  1:  Steady  Flow  around  a  Spur-Dyke 


The  developed  model  has  been  tested  using  the  measured  data 
of  Rajaratnam  and  Nwachukwu  (1983)  for  the  steady  flow 
around  a  spur-dyke  in  a  straight  tilting  rectangular  flume.  The 
flume  was  37  m  long  and  0.92  m  wide.  The  experimental  run  Al 
is  simulated  here,  in  which  the  flume  bed  and  walls  were 
smooth,  the  spur-dyke  was  a  thin  aluminum  plate  of  0.152  m  in 
length,  the  flow  discharge  was  0.0453  mVs,  and  the  approach 
flow  depth  was  0.189  m. 

To  check  whether  the  model  results  depend  on  the 
computational  mesh,  simulations  are  carried  out  using  two 
quadtree  meshes  and  a  uniform  mesh.  The  two  quadtree  meshes 
around  the  spur-dyke  are  shown  in  Figure  3,  and  both  consist  of 
three-level  refinements  near  the  spur-dyke.  The  finest  grids 
cover  the  entire  recirculation  zone  in  quadtree  mesh  A,  but  only 
the  region  near  the  spur-dyke  in  quadtree  mesh  B.  The  finest 
grid  spacing  near  the  spur-dyke  in  both  quadtree  meshes  is  about 
0.01  m  in  x  andy  directions.  For  comparison,  a  uniform  mesh 
with  this  fine  resolution  covering  the  entire  computational 
domain  is  also  tested.  The  number  of  cells  in  the  uniform  mesh, 
quadtree  mesh  A  and  B  is  57120,  19335,  and  6482,  respectively. 
All  model  parameters  are  set  the  same  for  the  three  simulation 
runs.  The  modified  mixing-length  turbulence  model  is  used  here, 
with  the  coefficient  cm  calibrated  as  0.3.  The  Manning’s  n  is  set 
as  0.012. 

Figure  4  shows  the  flow  patterns  calculated  using  the 
developed  model  with  the  three  meshes.  For  better  view,  certain 
points  are  skipped  on  the  plots.  One  can  see  that  the 
recirculation  and  main  flow  patterns  simulated  using  the  three 
meshes  are  quite  similar.  The  recirculation  zone  simulated  using 
quadtree  mesh  B  is  slightly  shorter  than  those  using  the  uniform 
mesh  and  quadtree  mesh  A.  Figure  5  compares  the  measured 
and  calculated  depth-averaged  flow  velocities  in  cross-sections 
located  at  x/b= 2,  4,  6  and  8.  Here,  b  is  the  length  of  the  spur- 
dyke.  The  model  reasonably  predicts  the  main  and  recirculation 


Figure  3.  Meshes  near  a  spur-dyke  in  a  straight  flume:  (a)  Quadtree 
mesh  A  and  (b)  Quadtree  mesh  B  (Circles:  locations  of  cell  centers). 


Figure  4.  Calculated  flow  patterns  near  the  spur-dyke  using  (a)  Uniform 
mesh,  (b)  Quadtree  mesh  A,  and  (c)  Quadtree  mesh  B. 
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Figure  5.  Comparison  of  measured  and  simulated  flow  velocities  at 
selected  cross-sections. 


Figures  6  and  7  show  the  contours  of  velocity  magnitude  and 
water  surface  elevation  simulated  using  quadtree  mesh  B.  The 
recirculation  flows  cover  several  levels  of  refined  mesh,  and  the 
transition  between  flows  on  fine  and  coarse  grids  is  very 
smooth.  This  demonstrates  that  the  numerical  discretization  and 
solution  methods  in  the  developed  model  are  adequate  to  handle 
the  unstructured  quadtree  mesh. 
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Figure  6.  Calculated  velocity  contours. 


Case  2:  Tidal  Flow  in  Gironde  Estuary 

The  Gironde  Estuary  is  located  in  southwestern  France.  It 
receives  runoff  from  the  Garonne  River  and  the  Dordogne  River 
and  empties  into  the  Atlantic  Ocean  (Figure  8).  The  water- 
surface  width  varies  from  2  to  14  km,  and  the  flow  depth  in  the 
navigation  channels  is  about  6-30  m.  The  estuary  is  partially 


x(m) 


Figure  7.  Calculated  water  level  contours. 


mixed  and  macrotidal,  with  a  12  hour  and  25  minutes  tidal  lunar 
period  and  a  tidal  amplitude  of  1.5-5  m  at  the  mouth  (Li  et  al., 
1994).  The  simulation  domain  is  about  80  km  long,  starting  from 
the  mouth  to  the  Garonne  River  and  the  Dordogne  River. 

The  bed  topography  was  provided  on  a  uniform  mesh,  with  a 
size  of  250  m  x  125  m  for  each  cell.  Considering  the  domain  is 
relatively  simple,  the  uniform  mesh  for  the  topography  is  used 
here  as  the  computational  mesh.  It  is  treated  as  the  simplest 
quadtree  mesh,  so  that  the  developed  model  can  handle  it 
straightforwardly.  The  data  measured  from  May  19  to  25,  1975 
is  used  to  validate  the  model.  The  computational  time  step  is  30 
minutes.  At  the  estuary  mouth,  the  tidal  elevation  is  given  by  the 
recorded  time  series  at  the  station  “Pointe  de  Grave”.  At  the  two 
upstream  ends,  the  flow  discharges  of  the  Garonne  River  and  the 
Dordogne  River  are  specified  according  to  the  measured  data  at 
La  Reole  and  Pessac.  The  Manning’s  n  is  set  as  0.015.  The 
modified  mixing-length  model  is  used  here,  with  the  coefficient 
cm  set  as  1 .2. 


Figure  8.  Sketch  of  Gironde  Estuary,  France. 
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The  flow  fields  in  flood  and  ebb  tides  are  reasonably  well 
predicted,  as  shown  in  Figure  9.  Figure  10  compares  the 
measured  and  simulated  water  levels  at  stations  lie  Verte  and 
Richard.  The  amplitude  and  phase  are  well  predicted  by  the 
numerical  model.  No  obvious  phase  difference  exists  between 
the  measured  and  simulated  tidal  levels.  Figure  11  shows  the 
comparison  of  the  measured  and  simulated  flow  velocities  at 
stations  Pauillac  and  Richard.  The  measured  flow  velocities  are 
1  m  under  the  water  surface  and  1  m  above  the  river  bed, 
respectively.  The  simulated  depth-averaged  flow  velocity  stays 
between  them.  The  agreement  is  reasonably  good. 


Figure  9.  Simulated  flow  patterns  in  Gironde  Estuary. 


Figure  10.  Measured  and  simulated  water  levels  at  selected  stations  in 
Gironde  Estuary. 


Figure  11.  Measured  and  simulated  velocities  at  selected  stations  in 
Gironde  Estuary. 


Case  3:  Tidal  Flow  in  Grays  Harbor 

Grays  Harbor  is  located  on  the  southwest  Washington  coast, 
USA,  at  the  mouth  of  the  Chehalis  River.  The  estuary  is  one  of 
the  largest  in  the  continental  United  States.  As  part  of  the  U.  S. 
Anny  Corps  of  Engineers  Grays  Harbor  Estuary  Physical 
Dynamics  Study,  current  and  waves  were  measured  during 
September  to  November  of  1999  (Osborne  et  al.,  2003). 
Figure  12  shows  the  plan  view  of  the  harbor  and  the  deployed 
measurement  stations.  Figure  13  shows  part  of  the  quadtree 
mesh  near  the  entrance  of  Grays  Harbor.  The  mesh  is  refined 
around  the  jetties  and  near  the  channels,  and  consists  of  96,100 
cells.  The  finest  grid  spacing  is  25  m  near  the  jetties  and  the 
coarsest  one  is  800  m  near  the  offshore  boundary  at  deep  water. 
The  computational  time  step  is  30  min.  The  measured  water 
levels  from  the  station  nearest  to  the  offshore  boundary  are  used 
as  the  boundary  condition.  The  wave  radiation  stresses 
calculated  by  the  CMS-Wave  model  are  considered  in  the 
simulation  of  current. 
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Figure  12.  Measurement  stations  deployed  in  Grays  Flarbor  in 
September  1999. 
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Figure  13.  Computational  mesh  near  the  entrance  of  Grays  Flarbor  (Dots 
indicate  cell  center  locations  and  lines  are  bed  elevation  contours). 
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Figure  14  shows  the  simulated  streamlines  during  ebb  tide 
near  the  entrance  of  Grays  Harbor.  One  can  see  the  streamlines 
smoothly  pass  through  coarse  and  fine  mesh  zones.  Figure  15 
shows  the  simulated  wetting  and  drying  processes  on  the 
floodplain  in  the  northeast  corner  of  the  bay  due  to  tide  change. 
Part  of  the  nodes  are  wetted  during  flood  tide  but  become  dry 
during  ebb  tide.  This  demonstrates  that  the  model  handles 
wetting  and  drying  efficiently.  Figure  16  compares  the 
computed  and  measured  water  levels  at  Stations  4  and  5,  and 
Figure  17  compares  the  computed  and  measured  current  speeds 
at  Stations  2  and  3  for  a  period  of  six  days  in  the  late  September 
of  1999.  The  model  has  reproduced  the  tidal  water  level  and 
current  speed  reasonably  well. 
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Figure  14.  Simulated  streamlines  during  ebb  tide  near  the  entrance  of 
Grays  Harbor  (Background  is  bed  elevation  contours). 


CONCLUSIONS 

An  implicit  finite  volume  scheme  with  several  advanced 
techniques  well  proven  in  CFD  has  been  applied  to  solve  the 
depth-averaged  2-D  shallow  water  flow  equations  to  improve 
the  computational  efficiency  of  the  existing  CMS  flow  model. 
This  new  model  uses  an  unstructured  multiple-level  quadtree 
(telescoping)  rectangular  mesh,  which  can  locally  refine  the 
mesh  around  structures  or  in  high-gradient  regions  and  thus 
improve  the  accuracy  of  the  model  with  a  relatively  small 
increase  in  number  of  cells.  The  grid  points  are  numbered  by 
means  of  an  unstructured  index  system,  so  that  the  mesh  can  be 
flexible  while  it  has  the  merits  of  the  traditional  rectangular 
mesh.  The  model  uses  the  non-staggered  (collocated)  system,  in 
which  primary  variables  v-velocity,  and  water  level  are 
stored  on  the  same  set  of  grid  points,  so  that  the  interface 
between  fine  and  coarse  cells  can  be  easily  handled. 

The  convective  tenns  in  the  momentum  equations  are 
discretized  using  one  of  several  upwinding  schemes,  including 
the  hybrid  upwind/central  difference,  exponential  difference  and 
HLPA  schemes.  The  HLPA  scheme  has  second  order  accuracy 
in  space,  while  the  other  two  schemes  have  the  accuracy 
between  first  and  second  orders.  The  SIMPLEC  algorithm  with 
under-relaxation  is  used  to  handle  the  coupling  of  water  level 
and  velocity  and  achieve  high  numerical  stability  and  efficiency. 
Fluxes  at  cell  faces  are  determined  using  the  Rhie  and  Chow’s 


Figure  15.  Floodplain  wetting  and  drying  due  to  tide. 


momentum  interpolation  method,  which  eliminates  spurious 
oscillations  existing  on  non-staggered  grid  that  can  occur  with 
linear  interpolation. 

Because  the  mesh  is  unstructured,  the  discretized  algebraic 
equations  have  sparse  matrices  of  coefficients.  These  equations 
are  solved  iteratively  using  the  flexible  GMRES  method  with 
ILUT  preconditioning,  which  is  efficient  for  solving  such 
algebraic  equations  (Saad,  1993). 

The  first  test  case  presented  in  this  paper  has  demonstrated 
that  the  developed  numerical  algorithms  handle  quadtree  meshes 
well  and  the  simulated  flow  structures  transit  continuously 
between  coarse  and  fine  cells.  The  second  and  third  test  cases 
were  used  to  test  the  developed  model’s  robustness  and 
reliability,  showing  that  a  time  step  of  about  a  half  hour  can  be 
used  in  simulation  of  tidal  flow  with  a  12-hr  period  and  the 
simulated  tidal  levels  and  velocities  are  in  good  agreement  with 
the  measured  data.  The  model  is  able  to  handle  the  wetting  and 
drying  problem  well  and  has  the  potential  to  be  used  in 
simulation  of  long-term  processes  near  coastal  inlets. 
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Figure  16.  Measured  vs.  computed  water  surface  elevations  at  Stations  4  and  5. 


Figure  17.  Measured  vs.  computed  current  speeds  at  Stations  2  and  3. 
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APPENDIX:  NOTATIONS 

am,  ctEk,  <*sk,  <*Nk,  tip  =  coefficients  in  discretized  equations 
a"  ,  a'k ,  ak  =  coefficients  for  u-,  v-  and  f  equations 

cm  =  empirical  coefficient  for  eddy  viscosity 

fc  =  Coriolis  force  coefficient 

fx.p’fv.p  =  weighting  factors  for  interpolation 

Fek,  F„k,  Fnk,  Fsk=  fluxes  at  cell  faces 

g  =  gravitational  acceleration 

h  =  total  flow  depth 

n  =  Manning’s  roughness  coefficient 

S  =  source  tenn 

t  =  time 

u,  v  =  depth-averaged  flow  velocities  in  x-  andy-directions 

u,  —  bed  shear  velocity 

x, y  =  horizontal  coordinates 
a„ ,  o.v  =  relaxation  factors 

AAP=  area  of  the  control  volume  at  node  P 

At  =  time  step  length 

Ay,  Ay  =  lengths  of  cell  face 

/;  =  water  level 

f  =  water  level  correction 

v,  =  eddy  viscosity 
p  =  flow  density 

Tbx,  tby  =  bed  shear  stresses 

subscript  k  =  index  of  the  faces  at  each  side 

subscript  /  =  index  of  neighboring  cells 

subscripts  w,  e,s,n  =  west,  east,  south,  north  sides  of  cell 

subscripts  W,  E,  S,  N=  west,  east,  south,  north  adjacent  nodes 

superscript  n  =  time  step  index 

superscript  o  =  value  of  previous  iteration  step 

superscript  *  =  guessed  or  approximate  value 


Journal  of  Coastal  Research,  Special  Issue  No.  59,  2011 


